Map Data

Global Administrative Areas 에서 다운로드한 zip 파일을 풀고 shape 파일을 mapshaper 에서 간략히 한 후 작업에 착수한다. .shp 형태로만 불러들일 수 있기 때문에 @data는 원본에서 복사하여야 한다.

Simplified Versions (from mapshaper)

시도

library(sf)
library(rgeos)
library(maptools) #> `readShapePoly()`, `rgdal::readOGR()`로 대체
library(GISTools) #> `choropleth`
library(ggmap) #> `geocode()`, `ggmap()`, `qmap()`, `revgeocode()`
library(ggplot2) #> `ggplot()`
library(rgdal) #> `CRS`, `ogrInfo()`, `ogrListLayers()`, `readOGR()`, `spTransform()`
library(dplyr) #> `arrange()`, `filter()`, `group_by()`, `left_join()`, `mutate()`,  
library(broom) #> `tidy()`
library(extrafont) #> "HCR Dotum LVT"
options(width = 132)

Simplified Shape files

shp_S <- readOGR(dsn = "../data/KOR_adm/", 
                layer = "KOR_adm1",
                stringsAsFactors = FALSE)
shp_N <- readOGR(dsn = "../data/PRK_adm/", 
                layer = "PRK_adm1",
                stringsAsFactors = FALSE)
ogrInfo(dsn = "../data/KOR_adm1",  #> `ogrInfo()`
        layer = "KOR_adm1")
shp_S_1 <- readOGR(dsn = "../data/KOR_adm1/", 
                layer = "KOR_adm1",
                stringsAsFactors = FALSE)
#                 encoding = "CP949")
ogrInfo(dsn = "../data/PRK_adm1",  #> `ogrInfo()`
        layer = "PRK_adm1")
shp_N_1 <- readOGR(dsn = "../data/PRK_adm/", 
                layer = "PRK_adm1",
                stringsAsFactors = FALSE)
#                 encoding = "CP949")

shp_S_1@data <- shp_S@data
shp_N_1@data <- shp_N@data
data.frame(shp_S_1)
data.frame(shp_N_1)
class(shp_S_1)
class(shp_N_1)
summary(shp_S_1)
summary(shp_N_1)
coordinates(shp_S_1)
coordinates(shp_N_1)
Names_S <- shp_S_1$NL_NAME_1
Names_N <- shp_N_1$NL_NAME_1
shp_S_1$Names_kr <- c("부산", "충북", "충남", "대구", "대전", "강원", "광주", "경기", "경북",
                      "경남", "인천", "제주", "전북", "전남", "세종", "서울", "울산")
shp_N_1$Names_kr <- c("자강", "함북", "함남", "황북", "황남", "개성", "강원", "금강", "평북",
                      "평남", "평양", "라선", "량강", "신의")
shp_S_1$Names_kr_old <- c("경상", "충청", "충청", "경상", "충청", "강원", "전라", "경기",
                          "경상", "경상", "경기", "전라", "전라", "전라", "충청", "서울",
                          "경상")
shp_N_1$Names_kr_old <- c("평안", "함길", "함길", "황해", "황해", "유후", "강원", "강원",
                          "평안", "평안", "평안", "함길", "함길", "평안")
shp_S_1$rates <- c(98.9, 33.3, 33.3, 98.9, 33.3, 12.0, 99.1, 98.6, 98.9, 98.9, 98.6, 
                   99.1, 99.1, 99.1, 33.3, 50.7, 98.9)
shp_N_1$rates <- c(4.5, 1.0, 1.0, 22.3, 22.3, 94.1, 12.0, 12.0, 4.5, 
                   4.5, 4.5, 1.0, 1.0, 4.5)
index_S <- shp_S_1$Names_kr %in% c("경북", "대전", "강원", "전북", "경기", "서울")
index_N <- shp_N_1$Names_kr %in% c("개성", "황북", "평남", "함남")
Lon_S <- coordinates(shp_S_1)[, 1]
Lat_S <- coordinates(shp_S_1)[, 2]
Lon_S_1 <- Lon_S[index_S]
Lat_S_1 <- Lat_S[index_S]
Lat_S_1[c(3, 6)] <- Lat_S_1[c(3, 6)] + c(-0.3, 0.2)
Lon_N <- coordinates(shp_N_1)[, 1]
Lat_N <- coordinates(shp_N_1)[, 2]
Lon_N_1 <- Lon_N[index_N]
Lat_N_1 <- Lat_N[index_N]
labels_S <- paste(shp_S_1$Names_kr_old[index_S], 
                  paste(shp_S_1$rates[index_S], "%", sep = ""), sep = ":")
labels_N <- paste(shp_N_1$Names_kr_old[index_N], 
                  paste(shp_N_1$rates[index_N], "%", sep = ""), sep = ":")
col_S <- c("black", "black", "white", "white", "white", "white")
col_N <- c("black", "black", "white", "black")

Simple Plots

library(extrafont)
par(mar = c(0, 0, 0, 0))
plot(shp_S_1, ylim = c(32, 44))
plot(shp_N_1, add = TRUE)
p_S <- pointLabel(Lon_S, Lat_S, Names_S, offset = 0, cex = 0.5, family = "HCR Dotum LVT")
p_N <- pointLabel(Lon_N, Lat_N, Names_N, offset = 0, cex = 0.5, family = "HCR Dotum LVT")

# png("./pics/Sejong_map.png")
# shades <- shading(breaks = seq(0, 100, by = 0.1), 
#                   cols = colorRampPalette(brewer.pal(11, "Spectral"))(1001))
shades <- shading(breaks = seq(0, 100, by = 0.1), 
                  cols = colorRampPalette(brewer.pal(11, "RdYlBu"))(1001))
# shades <- shading(breaks = seq(0, 100, by = 0.1), 
#                  cols = colorRampPalette(c("Red", "Pink", "Green", "Cyan"))(1001))
# rates_S_shades <- auto.shading(shp_S_1$rates, cols = rev(brewer.pal(5, "Greens")))
# rates_N_shades <- auto.shading(shp_N_1$rates, cols = rev(brewer.pal(5, "Greens")))
# choropleth(shp_S_1, shp_S_1$rates, shading = rates_S_shades, ylim = c(32, 44))
# choropleth(shp_N_1, shp_N_1$rates, shading = rates_N_shades, add = TRUE)
par(mar = c(0, 0, 0, 0), family = "HCR Dotum LVT")
choropleth(shp_S_1, shp_S_1$rates, shading = shades, ylim = c(32, 44))
choropleth(shp_N_1, shp_N_1$rates, shading = shades, add = TRUE)
text(Lon_S_1, Lat_S_1, labels = labels_S, col = col_S, cex = 0.8)
text(Lon_N_1, Lat_N_1, labels = labels_N, col = col_N, cex = 0.8)

# choro.legend(1110000, 1600000, AREA.shades)
# dev.off()
# dev.copy(png, file = "../pics/Sejong_ref_map.png", width = 960, height = 960)
# dev.off()

ggplot

shp_S_df <- tidy(shp_S_1)
## Regions defined for each Polygons
shp_N_df <- tidy(shp_N_1)
## Regions defined for each Polygons
g0 <- ggplot() +
  geom_polygon(data = shp_S_df, aes(x = long, y = lat, group = group),
               color = "black", fill = "white") +
  geom_polygon(data = shp_N_df, aes(x = long, y = lat, group = group),
               color = "black", fill = "white") 
g0 + 
  coord_quickmap()

Piece by piece into Unification

index_SL <- shp_S_1$Names_kr_old == "서울"
shp_SL <- shp_S_1[index_SL, ]
plot(shp_SL)

index_YH <- shp_N_1$Names_kr_old == "유후"
shp_YH <- shp_N_1[index_YH, ]
plot(shp_YH)

index_GS <- shp_S_1$Names_kr_old == "경상"
index_GS
##  [1]  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
shp_GS <- shp_S_1[index_GS, ]
shp_GS_u <- gUnaryUnion(shp_GS)
plot(shp_GS_u)

index_JL <- shp_S_1$Names_kr_old == "전라"
index_JL
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE FALSE
shp_JL <- shp_S_1[index_JL, ]
shp_JL_u <- gUnaryUnion(shp_JL)
plot(shp_JL_u)

index_CC <- shp_S_1$Names_kr_old == "충청"
index_CC
##  [1] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
shp_CC <- shp_S_1[index_CC, ]
shp_CC_u <- gUnaryUnion(shp_CC)
plot(shp_CC_u)

index_GG <- shp_S_1$Names_kr_old == "경기"
index_GG
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
shp_GG <- shp_S_1[index_GG, ]
#> Using a zero-width buffer cleans up many toplogical problems in R 
shp_GG <- gBuffer(shp_GG, byid = TRUE, width = 0)
## Warning in gBuffer(shp_GG, byid = TRUE, width = 0): Spatial object is not projected; GEOS expects planar coordinates
plot(shp_GG)

shp_GG_u <- gUnaryUnion(shp_GG)
plot(shp_GG_u)

index_PA <- shp_N_1$Names_kr_old == "평안"
index_PA
##  [1]  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE FALSE FALSE  TRUE
shp_PA <- shp_N_1[index_PA, ]
shp_PA_u <- gUnaryUnion(shp_PA)
plot(shp_PA_u)

index_GW <- shp_N_1$Names_kr_old == "강원"
index_GW
##  [1] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
shp_GW <- shp_N_1[index_GW, ]
shp_GW_u <- gUnaryUnion(shp_GW)
plot(shp_GW_u)

index_HH <- shp_N_1$Names_kr_old == "황해"
index_HH
##  [1] FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
shp_HH <- shp_N_1[index_HH, ]
shp_HH_u <- gUnaryUnion(shp_HH)
plot(shp_HH_u)

index_HG <- shp_N_1$Names_kr_old == "함길"
index_HG
##  [1] FALSE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE
shp_HG <- shp_N_1[index_HG, ]
shp_HG_u <- gUnaryUnion(shp_HG)
plot(shp_HG_u)

shp_GW_S <- shp_S_1[shp_S_1$Names_kr_old == "강원", ]
plot(shp_GW_u)
plot(shp_GW_S, add = TRUE)

shp_GW_I <- gIntersection(shp_GW_u, shp_GW_S)
class(shp_GW_I)
## [1] "SpatialCollections"
## attr(,"package")
## [1] "rgeos"
plot(shp_GW_I)

# shp_GW_u_sf <- st_as_sfc(shp_GW_u)
# class(shp_GW_u_sf)
# plot(shp_GW_u_sf)
# shp_GW_S_sf <- st_as_sfc(shp_GW_S)
# class(shp_GW_S_sf)
# plot(shp_GW_S_sf)
# shp_GW_st_u <- st_union(shp_GW_u_sf, shp_GW_S_sf)
# plot(shp_GW_st_u)
# shp_GW_st_i <- st_intersection(shp_GW_u_sf, shp_GW_S_sf)
# plot(shp_GW_st_i)
# shp_GW_st_d <- st_sym_difference(shp_GW_st_u, shp_GW_st_i)
# plot(shp_GW_st_d)
# shp_GW_st_c <- st_union(shp_GW_st_d)
# plot(shp_GW_st_c)
shp_GW <- gUnion(shp_GW_u, shp_GW_S)
class(shp_GW)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
coordinates(shp_GW)
##       [,1]     [,2]
## 1 127.9615 38.10755
coordinates(shp_GW_u)
##       [,1]     [,2]
## 1 127.4634 38.71698
coordinates(shp_GW_S)
##       [,1]    [,2]
## 5 128.2809 37.7168
# shp_GW_spdf <- SpatialPolygonsDataFrame(shp_GW, 
#                                         data = data.frame(dummy = 1:nrow(coordinates(shp_GW))))
# writeOGR(shp_GW_spdf, dsn = "../data/GW", layer = "GW", driver = "ESRI Shapefile")
plot(shp_GW)

summary(shp_GW)
## Object of class SpatialPolygons
## Coordinates:
##         min       max
## x 126.67380 129.35239
## y  37.00888  39.42428
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
class(shp_GW_u)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"
class(shp_CC)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

Old Map

par(mar = c(0, 0, 0, 0))
plot(shp_HG_u, xlim = c(123, 132), ylim = c(32, 44))
plot(shp_PA_u, add = TRUE)
plot(shp_HH_u, add = TRUE)
plot(shp_GW, add = TRUE)
plot(shp_GG_u, add = TRUE)
plot(shp_CC_u, add = TRUE)
plot(shp_JL_u, add = TRUE)
plot(shp_GS_u, add = TRUE)
plot(shp_SL, add = TRUE)
plot(shp_YH, add = TRUE)

Colours

cols <- unique(c(shp_S_1$rates, shp_N_1$rates)) * 10
pie(rep(1, 10), col = colorRampPalette(brewer.pal(11, "RdYlBu"))(1001)[cols])

names(cols) <- c("GS", "CC", "GW", "JL", "GG", "SL", "PA", "HG", "HH", "YH")
cols_hex <- colorRampPalette(brewer.pal(11, "RdYlBu"))(1001)[cols]
names(cols_hex) <- names(cols)
cols_hex
##        GS        CC        GW        JL        GG        SL        PA        HG        HH        YH 
## "#333D98" "#FDBE70" "#DC3B2C" "#333C98" "#343F99" "#FDFEC2" "#BB1526" "#A90426" "#F57B49" "#3C5BA7"

Combined

par(mar = c(0, 0, 0, 0))
plot(shp_HG_u, xlim = c(123, 132), ylim = c(32, 44), col = cols_hex["HG"])
plot(shp_PA_u, add = TRUE, col = cols_hex["PA"])
plot(shp_HH_u, add = TRUE, col = cols_hex["HH"])
plot(shp_GW, add = TRUE, col = cols_hex["GW"])
plot(shp_GG_u, add = TRUE, col = cols_hex["GG"])
plot(shp_CC_u, add = TRUE, col = cols_hex["CC"])
plot(shp_JL_u, add = TRUE, col = cols_hex["JL"])
plot(shp_GS_u, add = TRUE, col = cols_hex["GS"])
plot(shp_SL, add = TRUE, col = cols_hex["SL"])
plot(shp_YH, add = TRUE, col = cols_hex["YH"])